Study on the freeze–thaw effect on concrete arch dam using an improved response surface method

The damming materials will be deteriorated during the long-term service of concrete arch dam by weathering, seepage dissolution, freeze–thaw and other factors which will put dams at potential risk. However, there are few researches on the reliability assessment of arch dams considering the effect of material deterioration. In this paper, a prediction model of concrete mechanical properties under freeze–thaw cycles is established using regression analysis. Considering various uncertain factors, a time-varying structural reliability analysis model is established using finite element equivalent stress method and an optimized response surface method with improved sampling strategy and convergence criterion. A case study is later carried out. The results show that, at the early stage of freeze–thaw, the structural failure probability will increase steadily in a small range. However, as the tensile strength of concrete continues to decrease under long-time freeze–thaw action, the reliability of the arch dam structure will suddenly fall off a cliff. The proposed deterioration model and reliability analysis probability model are reasonable and practical and could be serve as a technical support for similar engineering projects.

mechanical uncertainty and so on.The structural reliability research is relatively scarce.Yilmaz et al. proposed a stochastic deterioration modelling approach to evaluate the performance deterioration of corroded concrete structures caused by reinforcement corrosion 15 .However, it does not reflect the effect of concrete deterioration on arch dam.In general, the reliability of the arch dam structure under freeze-thaw effect, considering various uncertainty factors, still needs further study.
In this paper, a freezing-thawing deterioration material model of concrete is proposed through regression analysis.The time-varying strength and reliability of an arch dam under freezing-thawing action are analyzed by using a finite element equivalent stress method and response surface method.The research results provide a method and basis for the safety evaluation of arch dam structure considering material deterioration.Other material deterioration progress, such as weathering, seepage dissolution, et al., could also be similarly analyzed by the proposed method.

Deterioration model of concrete under freeze-thaw effect
freeze-thaw action is generally considered to be the internal microstructure changing process of cement mortar due to repeated freezing and melting of internal moisture.The freeze-thaw action mainly affects the strength and elastic modulus of concrete.Under freeze-thaw effect, the dynamic elastic modulus of concrete E D gradually decreases with the increase of freeze-thaw times, and the change rule is basically the same.In order to study the time-varying reliability process of arch dam under static condition, the dynamic elastic modulus E D needs to be transformed into static elastic modulus E S .According to "Standard for seismic design of hydraulic structures" (GB 51247-2018) 16 , under seismic conditions, the dynamic elastic modulus of concrete increases by 50% on the basis of static elastic modulus.It could be assumed that the dynamic elastic modulus of concrete under freeze-thaw action remains the same as 1.5 times of static elastic modulus, namely: Based on the exponential fitting model regression analysis of the dynamic elastic modulus test data in relevant literatures 2,[17][18][19][20][21][22][23] , the test data of plain concrete are extracted.And the static elastic modulus of concrete under freeze-thaw action could be reverse analyzed and obatined which is shown in Fig. 1: where E S,0 , E D,0 are the initial static and dynamic elastic modulus of concrete without freeze-thaw action, respectively; E S,n , E D,n are the static and dynamic elastic modulus of concrete after n times of freeze-thaw action.
Based on a set of freeze-thaw cycle tests, Liu 22 also proposed a deterioration model which is shown in Eq. (3) and Fig. 1.
It could be noted that within the first 75 freeze-thaw cycles, the results of the model proposed in this paper are not much different from Liu's model.However, the failure of Liu's concrete strength test after 200 cycles caused the regression model to decline rapidly after 75 cycles.As more experimental data are used in this paper, a smoother deterioration characteristic curve is obtained.
According to research results and test data in relevant literature [23][24][25][26] , the relative residual compressive strength of concrete under freeze-thaw action has a good correlation with the relative dynamic elastic modulus and follows the power function relationship.By fitting the data in the literature through power funtion, the relation between the residual compressive strength of concrete and the relative dynamic elastic modulus could be determined and shown in Fig. 2: (1) where f c,0 is the initial compressive strength of concrete without freeze-thaw action; f c,n is the compressive strength of concrete after n times of freeze-thaw action.
It could be noted that the deterioration curve proposed in this paper lies in the middle of other models.Zhao's test data did not include cases where the relative residual dynamic elastic modulus is below 0.9 which makes the regressed deterioration curve much steeper when the relative residual dynamic elastic modulus decreases further.Hong thinks when damage degree is less than 40%, the compressive strength of freeze-thaw damaged concrete is not lower than 70% of its original strength.Therefore, the conservative estimate of the residual compressive strength of freeze-thaw damaged concrete is equivalent to the residual dynamic elastic modulus.As more experimental data are used in this paper, a relative compromised compressive strength deterioration characteristic curve is obtained.
Assuming that the proportional relationship between compressive strength and tensile strength of concrete under freeze-thaw action is the same as the standard value of concrete strength 27 , the tensile strength of concrete under freeze-thaw action could be determined as follows: where f c,n , f t,n are the compressive strength and tensile strength of concrete after n times of freeze-thaw action, respectively; f ck , f tk are the standard compressive strength and tensile strength of concrete, respectively.
Equations ( 2)-( 4) define respectively the static elastic modulus and strength deterioration process of concrete under freeze-thaw action.

Finite element equivalent stress method
An input-output model is needed to evaluate the reliability of arch dam structure.In the input-output model defined by finite element method, stress concentration could occur at the corner edge of the dam resulting in the decrease of the reliability index 28 .In order to reduce stress concentration, the finite element internal force analysis results will be treated by equivalent stress method in the post-processing stage 29 .
The stress at each node of arch dam σ ij ′ in the global coordinate system (x′, y′, z′) is first calculated by finite element method.Then, a local coordinate system (x, y, z) with radial and tangential axes of the arch is established at node k of the horizontal arch circle, where the X-axis is parallel to the tangential direction of the arch center line, the Y-axis is parallel to the radius direction, and the Z-axis is the vertical direction.The coordinate's origin is on the arch center line.The stress σ ij ′ in the global coordinate system could be converted into the stress in the local coordinate system σ ij = (σ x , σ y , σ z , τ xy , τ yz , τ zx ) T by means of the coordinate transformation matrix: where [T] is the coordinate transformation matrix.where (l i , m i , n i ) is the direction cosine of coordinates x, y and z.
On the horizontal section of the beam, the unit width on the center line of the arch is taken.The width at point y could be 1 + y/r.r is the curvature radius of the arch center line.The internal force of the beam could be obtained by integrating the stress and its moment along the thickness direction: where T is the thickness of arch circle, y 0 is the centroid coordinates of the beam section.
The radial section width of arch ring per unit height is defined as 1.The internal force of an arch could be obtained by integrating the arch stress and its moment along the thickness direction: Since the shear stresses are paired, the vertical shear and torque of the arch could be determined without calculation.Assuming that σ z , σ x , τ zx are distributed linearly from upstream to downstream, the stress of dam body could be calculated by the method of material mechanics according to the internal forces of arch and beam.The remaining stress components need to be obtained by taking the microelements from the upstream and downstream edges and solving the equilibrium equations of the microelements.

Downstream surface stress
The vertical normal stress of the cantilever beam on the horizontal plane could be calculated as follows: The horizontal normal stress of the arch in the radial plane could be calculated as follows: www.nature.com/scientificreports/ The tangential horizontal shear force of the horizontal cantilever beam could be calculated as follows: The other three stress components could be calculated as follows: where φ d is the angle between the dam surface and the vertical direction in the radial vertical plane; η d is the angle between the dam surface and the tangent line of the arch center line in the horizontal plane.
The three stresses parallel to the dam surface could be calculated as follows: The equivalent principal stress of the downstream dam surface could be calculated as follows: where where t is the arch ring thickness; A b = A a = t; R u = r + t/2; R d = r − t/2; P d is the downstream dam surface pressure.

Upstream surface stress
The upstream dam surface stress calculation is the same as the downstream.The main calculation formulas are as follows: (20) where P u is the upstream dam surface pressure.

Response surface method
When the probability distribution of input variables is determined according to the test data or specification parameters, the response surface method could be adopted to analyze the structural reliability 30 .In response surface method, a fitted response surface is adopted to replace the complex functional equation, which is essentially an improvement of Monte Carlo method and could greatly improve the efficiency of reliability evaluation 31 .However, the calculation accuracy of response surface method depends on the fitting degree of the selected sample points near the limit state surface 32 .It is necessary to obtain sufficient samples near the limit state surface.But too much sampling will affect the computational efficiency, and has little effect on the accuracy improvement.Therefore, an improved sampling strategy and convergence criterion are proposed to achieve a balance between computational accuracy and efficiency.The quadratic polynomial is usually selected as the response surface function form.Set X = (x 1 , x 2 , …, x n ) as the random input variable.The response surface function without cross terms could be defined as follows: As the Latin hypercube sampling algorithm is selected for sampling, at least 2n + 1 samples will be needed to fit the response surface function containing n random variables.However, 2n + 1 samples may not be close enough to the limit state surface, so the fitting effect could not be determined.The optimal response surface could be determined by stepwise trial calculation, which is expressed as follows: The ith step trial is defined to fit k response surfaces.The number of response sampling points used to fit each response surface is respectively in The target probability of each response face is solved separately as P f,i , P f,i+1 , P f,i+2 , …, P f,i+k-1 .If the following formula is true: where max(P f ), min(P f ), E(P f ) are respectively the maximum value, minimum value and mean value of the target probability in the set of stepped response surfaces, ε is the allowable error.
Then, it is considered that the ith step trial calculation has achieved a suitable fitting effect near the limit state surface.The response surface whose target probability is closest to E(P f ) will be taken as the representative fitting response surface.Otherwise, i + 1 step trial calculation will be conducted, and n sample points will be newly added on the basis of the original sample points for trial test until Formula (41) is satisfied.
In this method, the number of sample points is gradually expanded, and the fitting effect is constantly selfevaluated so as to determine the best sample point number and the corresponding response surface.The essence is to propose a computational convergence criterion.

Project overview
A concrete arch dam located in the northwest of China is studied.The maximum height is 81.6 m.The top elevation is 401.60 m, the bottom elevation is 320.00 m, the top width is 5.0 m, and the bottom width of the crown cantilever is 15.0 m.The dam was constructed amd completed in 1973 and has been in service for 50 years.The shape parameters of the arch dam are shown in Table 1.(35) The normal water level of the reservoir is 398.00 m, and the corresponding downstream water level is 341.70 m.The check flood level is 399.80 m, and the corresponding downstream water level is 343.00 m.The dead water level is 351.00 m, and the corresponding downstream water level is 326.50 m.The elevation of silt deposited in front of the dam is 328.00 m.
The temperature change process of the dam site is shown in Fig. 3.The average annual temperature at the dam site is 15.9 °C with an annual temperature amplitude of 16 °C.

Finite element model
The bedrock is argillaceous limestone which is hard and solid.Thus, the linear elastic model and hexahedral eight-node entity units are adopted to simulate the dam body and foundation.The calculation model is discretized into 101,228 nodes and 92,805 cells.The dam foundation model extends 2.5 times the length from upstream to downstream, 1.5 times the width from left bank to right bank of the dam abutment, and 2 times the arch height downward from the bottom of the arch dam (coordinate system: the origin is located at the vertex of the upstream surface of the arch beam, the X-axis is along the water flow direction, the upstream direction is negative, the downstream direction is positive, the Y-axis is along the vertical direction, the upward is positive, the downward is negative, the Z-axis is from right bank to left bank, the left bank is positive, the right bank is negative).Rigid constraints are applied to the bottom surface of the foundation, and linkage constraints are applied on the lateral surface along the direction of water flow.The overall model mesh is shown in Fig. 4. The arch dam body mesh is shown in Fig. 5.

Working conditions, loads and parameters
According to the actual operation of the arch dam and "Design specification for concrete arch dams (SL282-2018)" 33 , the working conditions of the dam strength reliability analysis could be determined and shown in Table 2.  Temperature load is one of the main loads of arch dam 34 .Recent studies show that the arch sealing temperature and residual temperature are the main sources of temperature load uncertainty 35 .In this paper, the arch sealing temperatures are treated as random variables to characterize the uncertainty of temperature load.The arch sealing temperatures of eight layers of arch rings from arch bottom to arch top were separately considered as independent random variables.The coefficient of variation is set as 0.05.Temperature load is calculated according to "Design specification for concrete arch dams (SL282-2018)" 33 considering average body temperature T m and equivalent linearly-distributed temperature difference T d .The temperature load when the arch sealing temperature is 16.0 °C is shown in Table 3.The thermal linear expansion coefficient is taken as 8 × 10 −6 .
According to laboratory test, the elastic modulus of dam body varies from 16.8 to 22.6GPa.The standard deviation is 2.28GPa.Thus, the median value of 19GPa is adopted and variation coefficient is taken as 2.28/19 = 0.12.Other material property parameters are generally taken in the same way.The Initial tensile strength of concrete defined according to Eq. ( 4).
The dam body is made of C15 concrete.According to the laboratory test, the recommended initial elastic modulus is 19GPa.The Poisson's ratio is 0.20.The density is 2350 kg/m 3 .The initial tensile strength is 2.74 MPa.
The foundation is argillaceous limestone.According to the laboratory test, the recommended elastic modulus is 6 GPa.The Poisson's ratio is 0.25.The density is 2600 kg/m 3 .
The arch sealing temperature is usually taken as the temperature that is a little higher than the average annual temperature at the dam site, which is 15.9 °C.Thus, the arch sealing temperature is taken as 16.0 °C.
Random variables and their statistic parameters are shown in Table 4.
Because the compressive strength of arch dam is generally rich, and arch dams are mostly damaged by tension, the failure mode is selected as tensile failure in the structural reliability analysis.www.nature.com/scientificreports/

Results and discussions
The initial mean value of each parameter is adopted to simulate the initial dam stress state.As shown in Fig. 6 (the finite element stress is positive in tensile stress and negative in compressive stress), the maximum principal tensile stress is 3.06 MPa and the maximum principal compressive stress is 6.80 MPa at the basic load combination condition, both of which occur in the normal water level and temperature rise condition (Working Condition 2).As the occurrence frequency of working condition ( 2) is much higher than other working conditions, this working condition becomes the actual control condition of arch dam reliability.The nephogram of principal stress shows that the finite element simulation produces local stress concentration.Thus, equivalent stress treatment is adopted and the results are shown in Table 5.The stress concentration is alleviated after treatment.
The response surface method proposed in this paper is adopted to evaluate the structural reliability.According to the definition of aggregation number k of stepped response surfaces expressed above, the larger the aggregation number k is, the harsher the convergence condition will be.Generally, the calculation results will be more secure, meanwhile the corresponding amount of calculation will be increased.In fact, a good convergence could be achieved when k ≥ 3. Considering the calculation accuracy and efficiency comprehensively, the aggregation number of stepped response surfaces k is set to be 4, and the allowable error ε is set to be 0.5.In order to directly reflect the cumulative probability distribution of structural stress, the maximum principal tensile stress is taken as the objective function.The initial concrete tensile strength is taken as a constant to solve the structural failure probability at initial state (without freeze-thaw effect).The stepwise response surface convergence process is shown in Table 6.The failure probability reaches relative stability in step 5, and the calculation effect has been guaranteed with 121 samples.According to the response surface selection strategy described above, response surface 7 is selected.The cumulative probability curve of maximum principal tensile stress output by this response surface is shown in Fig. 7.
It is assumed that a freeze-thaw cycle is completed when the daily minimum temperature is below − 1 °C for over five consecutive days followed by five consecutive days above -1 °C.According to the statistics of the reservoir operating temperature shown in Fig. 3, 29 freeze-thawing times have occurred so far.According to the freeze-thaw deterioration model of concrete proposed in this paper, the time-varying process of the dam material properties and arch dam failure probability under freeze-thaw action are calculated, as shown in Figs. 8 and 9.
According to the calculation results, it could be known that: (1) Compared with the dead water level working condition, the maximum principal tensile stress cumulative probability curve in the condition with higher water level ((1), ( 2), ( 5 3) With the increase of freeze-thawing times, the failure probability of arch dam structure generally increases and the reliability decreases.At early stage of freeze-thaw action (approximately within 50 times of freezethaw cycles), though the tensile strength of concrete is decreasing, the failure probability will increase steadily in a small range, and the reliability of the arch dam structure will not decrease significantly.Because the stress level of dam body is also relieved when the elastic modulus is reduced.After that, as the tensile strength of concrete continue to decrease, the structural failure probability will suddenly increase.After 80 times of freeze-thawing cycles, the annual failure probability of the arch dam structure will reach about 1% and the arch dam structure will be unreliable.(4) Although the residual strength decreases slowly due to freeze-thaw effect, it has a great influence on the reliability of the structural strength.Because the residual tensile strength slowly approaches the maximum tensile stress of the structure.As the dam has freeze-thaw action 29 times during its 50 years of service, it is estimated that the dam will reach a 1% probability of failure in about 88 years.The problem of freeze-thaw damage should be paid attention to and corresponding pre-strengthening measures should be taken.

Conclusions
The material of concrete arch dam is easy to deteriorate under freeze-thaw action which affects its structural reliability and service life.In this paper, a freeze-thaw deterioration model of mechanical properties of concrete materials is established based on measured data.A response surface method with improved convergence criterion is proposed.The time-varying reliability process of arch dam structure under freeze-thaw action is studied.The main conclusions could be made as follows: (1) The relationships between freeze-thaw cycles, concrete dynamic elastic modulus, static elastic modulus and compressive and tensile strength are analyzed.A freeze-thaw deterioration model of concrete is established through regression analysis to simulate the deterioration process of mechanical properties which could provide a simple and practical analysis method for engineering application.(2) In view of the contradiction between the computation effect and computation amount of response surface method, a convergence criterion and strategy for sampling is proposed so as to determine the relatively optimal response surface.The case study results show that the response surface fitted after the fifth step trial is basically stable, and the fluctuation is not obvious.To some extent, the proposed convergence criterion could guarantee the effectiveness of response surface calculation, avoid excessive sampling quantity, and reduce the amount of calculation.It is feasible to use the convergence criterion proposed in this paper to determine the relative optimal response surface.(3) The simulation results of time varying structural reliability of arch dam show that, at early of freezethaw, though the tensile strength of concrete is decreasing, the reliability of the arch dam structure will not decrease significantly, because the stress level of dam body is also relieved when the elastic modulus is reduced.However, as the tensile strength continues to decrease under freeze-thaw action, the failure probability will increase rapidly.After 50 times of freeze-thawing cycles, the structural reliability will suddenly fall off a cliff.As the dam has freeze-thaw action 29 times during its 50 years of service, it is estimated that the dam will reach a 1% probability of failure in about 88 years.The problem of freeze-thaw damage should be paid attention to and corresponding pre-strengthening measures should be taken.Learning from the stress distribution, the stress level of the arch end of both sides is high, and the weak parts such as the arch end of both sides are suggested being monitored and maintained in daily operation.

Figure 1 .
Figure 1.Reversed dynamic elastic modulus of concrete under freeze-thaw action.

Figure 3 .
Figure 3. Temperature change process at the dam site.

Figure 7 .
Figure 7. Cumulative probability curve of maximum principal tensile stress under initial condition.

Figure 8 .Figure 9 .
Figure 8. Time-varying process of the dam material properties under freeze-thaw action (mean value).

Table 1 .
Shape parameters of the arch dam.

Table 2 .
Dam structural reliability calculation conditions.

Table 3 .
Temperature load with the arch sealing temperature of 16.0 °C.

Table 4 .
Random variables and statistic parameters.The water depth corresponds to the normal storage level, dead level and check flood level respectively.

Table 5 .
Finite element equivalent stress results of working condition (2), unit: MPa.The maximum stress locations are basically the same as the location of direct finite element calculation results.

Table 6 .
Failure probability convergence process through response surface method.C(*) refers to working condition (*).
)) has lower slope and longer span, which indicates that the possible stress fluctuation range of arch dam body is larger and the arch structure is in a relatively unstable state under high water pressure load.(2)The calculation results of normal water level & temperature rise condition (working condition (2)) are close to that of dead water level & temperature drop condition.However, since the former happens almost every year after impoundment, it is the actual control condition of the failure probability of arch dam.The other working conditions occur less frequently in one year and have little influence on the annual failure probability.Therefore, the failure probability of arch dam is approximately equal to the failure probability of normal water level & temperature rise condition (working condition (2)).(